
********
* CPAP *
********

* DME Claims
{
use ".\statadata\dme_total_matched", clear
preserve
keep if inlist(hcpcs_cd, "E0601", "E0470", "E0471", "E0472") & mod1=="RR" & year(clm_thru_dt)==2010
gen CPAP = 1 if hcpcs_cd=="E0601"
collapse (sum) total_amt = line_alowd_chrg_amt q = line_srvc_cnt, by(CPAP)
egen total_q = sum(q)
gen share = q / total_q
restore
keep if hcpcs_cd=="E0601" & mod1=="RR" // CPAP
gen yq = qofd(clm_thru_dt)
format yq %tq
collapse (count) n_claim = line_alowd_chrg_amt (sum) total_amt = line_alowd_chrg_amt q = line_srvc_cnt, by(bene_id yq)
save ".\temp\sample_CPAP.dta", replace
}


* MBSF: Enrollment info
{
keep bene_id
duplicates drop
merge 1:m bene_id using ".\temp\mbsf", keep(3) nogen
keep bene_id year geo_id cba_round enroll* dual* 
reshape long enroll_Q dual_Q, i(bene_id year geo_id cba_round) j(qtr)
gen yq = yq(year,qtr)
format yq %tq

keep if enroll_Q==3 
gen dual = 1 if dual_Q==3
replace dual = 0 if dual_Q==0
bysort bene_id: egen temp = count(enroll_Q)
gen enroll_28 = 1 if temp==28 
keep bene_id yq geo_id cba_round dual enroll_28

merge 1:m bene_id yq using ".\temp\sample_CPAP.dta", keep(3) nogen
save ".\temp\sample_CPAP.dta", replace
}


*** Supplier ***
{
use ".\statadata\dme_total_matched", clear
keep if hcpcs_cd=="E0601" & mod1=="RR" // CPAP
gen yq = qofd(clm_thru_dt)
format yq %tq
gen year = year(dofq(yq))
merge m:1 bene_id year using ".\temp\mbsf", keepusing(geo_id) keep(1 3) nogen

collapse (count) n_claim = line_alowd_chrg_amt (sum) total_amt = line_alowd_chrg_amt, by(prvdr_npi yq geo_id)
bysort geo_id yq: egen total_total_amt = sum(total_amt)
gen share = total_amt/total_total_amt
gen share2 = share^2
gen n = 1 if share>0.01

collapse (sum) N_NPI = n HHI_NPI = share2 , by(geo_id yq)
drop if geo_id==""
drop if yq==.
save ".\temp\sample_CPAP_supplier", replace
}


*** Regression Sample Construction ***
{
use ".\temp\sample_CPAP.dta", clear

bysort bene_id geo_id: egen temp = min(yq)
cap gen n_initiator = 1 if temp==yq
gen n_bene = 1 if q!=.

gen n_bene2 = 1 if enroll_28==1
gen n_initiator2 = n_initiator if enroll_28==1
gen q2 = q if enroll_28==1
gen total_amt2 = total_amt if enroll_28==1
gen n_claim2 = n_claim if enroll_28==1

collapse (count) n_bene (sum) n_bene2 n_initiator* q* total_amt* n_claim*, by(yq geo_id cba_round)
drop if yq==.
drop if geo_id==""
unique yq
unique geo_id
di 969*28
save ".\temp\sample_CPAP_reg", replace

keep geo_id cba_round
duplicates drop
expand 28
bysort geo_id: gen yq = yq(2008,4)+_n
format yq %tq
merge 1:1 geo_id yq using ".\temp\sample_CPAP_reg.dta", nogen
merge 1:1 geo_id yq using ".\temp\sample_CPAP_supplier.dta", nogen
merge m:1 geo_id yq using ".\temp\population", keep(1 3) nogen

replace q = 0 if q==.
gen qX = q/pop_FFS*1000
gen lg_qX = log(qX)

replace n_bene = 0 if n_bene==.
gen n_beneX = n_bene/pop_FFS*1000
gen lg_n_beneX = log(n_beneX)

replace n_initiator = 0 if n_initiator==.
gen NinitiatorX = n_initiator/pop_FFS*1000
gen lg_NinitiatorX = log(NinitiatorX)	

replace total_amt = 0 if total_amt==.
gen paymentX = total_amt/pop_FFS*1000
gen lg_paymentX = log(paymentX)

gen price = total_amt/q
gen lg_price = log(price)

replace N_NPI = 0 if N_NPI==.
gen N_NPIX = N_NPI/pop_FFS*1000
foreach v of varlist N_NPI* HHI* {
	replace `v' = . if q==0
}

gen dd = 0
replace dd = 1 if yq>=yq(2011,1) & yq<=yq(2013,4) & cba_round==1
replace dd = 1 if yq>=yq(2013,3) & yq<=yq(2016,2) & cba_round==2
replace dd = 1 if yq>=yq(2014,1) & yq<=yq(2016,4) & cba_round==1
save ".\temp\sample_CPAP_reg.dta", replace
}


*** Regression: DID Main effect (CPAP) ***
use ".\temp\sample_CPAP_reg.dta", clear
egen geo_idnum = group(geo_id)

* Main
{
local rp replace
foreach ver in price qX n_beneX NinitiatorX paymentX {
	cap n qui reghdfe lg_`ver' dd [aweight=pop_FFS_2010], nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
foreach ver in HHI_NPI N_NPI {
	cap n qui reghdfe `ver' dd [aweight=pop_FFS_2010], nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("Table2_DIDCPAP") sheetreplace	
restore
}

* By Wave
{
local rp replace
foreach ver in price qX {
	cap n qui reghdfe lg_`ver' dd [aweight=pop_FFS_2010] if cba_round==0 | cba_round==1, nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Wave, R1, Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
foreach ver in price qX {
	cap n qui reghdfe lg_`ver' dd [aweight=pop_FFS_2010] if (cba_round==0 | cba_round==1) & yq<yq(2014,1), nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Wave, R1_nRC, Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
foreach ver in price qX {
	cap n qui reghdfe lg_`ver' dd [aweight=pop_FFS_2010] if cba_round==0 | cba_round==2, nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Wave, R2, Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA8_DIDCPAP_Wave") sheetreplace	
restore

local rp replace
foreach ver in HHI_NPI N_NPI {
	cap n qui reghdfe `ver' dd [aweight=pop_FFS_2010] if cba_round==0 | cba_round==1, nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Wave, R1, Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
foreach ver in HHI_NPI N_NPI {
	cap n qui reghdfe `ver' dd [aweight=pop_FFS_2010] if (cba_round==0 | cba_round==1) & yq<yq(2014,1), nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Wave, R1_nRC, Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
foreach ver in HHI_NPI N_NPI {
	cap n qui reghdfe `ver' dd [aweight=pop_FFS_2010] if cba_round==0 | cba_round==2, nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Wave, R2, Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA9_DIDCPAP_Wave") sheetreplace	
restore
}

* Geographic spillover
{
use ".\raw\ZIP_CBA_CBSA.dta", clear
keep geo_id CBA_adjacent 
duplicates drop
drop if geo_id==""

merge 1:m geo_id using ".\temp\sample_CPAP_reg.dta", keep(2 3) nogen
egen geo_idnum = group(geo_id)

local rp replace
foreach ver in price qX {
	cap n qui reg lg_`ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if CBA_adjacent==., vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010] if e(sample)==1
		cap n outreg2 using $resdir\reg.dta, dta keep(dd) addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
	cap n qui reg lg_`ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if CBA_adjacent==.
	est store r1_`ver'
	local r1_`ver'_dd = _b[dd]
}
foreach ver in HHI_NPI N_NPI {
	cap n qui reg `ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if CBA_adjacent==., vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010] if e(sample)==1
		cap n outreg2 using $resdir\reg.dta, dta keep(dd) addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
	cap n qui reg `ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if CBA_adjacent==.
	est store r1_`ver'
	local r1_`ver'_dd = _b[dd]
}

foreach ver in price qX {
	cap n qui reg lg_`ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if CBA_adjacent!=. | cba_round!=0, vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010] if e(sample)==1
		cap n outreg2 using $resdir\reg.dta, dta keep(dd) addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
foreach ver in HHI_NPI N_NPI {
	cap n qui reg `ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if CBA_adjacent!=. | cba_round!=0, vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010] if e(sample)==1
		cap n outreg2 using $resdir\reg.dta, dta keep(dd) addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}

replace dd = 1 if yq>=yq(2011,1) & yq<=yq(2013,4) & CBA_adjacent==1
replace dd = 1 if yq>=yq(2013,3) & yq<=yq(2016,2) & CBA_adjacent==2
replace dd = 1 if yq>=yq(2014,1) & yq<=yq(2016,4) & CBA_adjacent==1

foreach ver in price qX {
	cap n qui reg lg_`ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if cba_round==0
	est store r3_`ver'
	local r3_`ver'_dd = _b[dd]
	local diff_`ver'_dd = `r3_`ver'_dd'-`r1_`ver'_dd'
	
	qui suest r1_`ver' r3_`ver', vce(cluster geo_id)
	test [r1_`ver'_mean]dd = [r3_`ver'_mean]dd
	local diff_`ver'_dd_p = r(p)
	local diff_`ver'_dd_p : di %9.0g `diff_`ver'_dd_p'
	
	cap n qui reg lg_`ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if cba_round==0, vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010] if e(sample)==1
		cap n outreg2 using $resdir\reg.dta, dta keep(dd) addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean), Panel C - Panel A, `diff_`ver'_dd', P value, `diff_`ver'_dd_p') `rp'
	local rp
}
foreach ver in HHI_NPI N_NPI {
	cap n qui reg `ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if cba_round==0
	est store r3_`ver'
	local r3_`ver'_dd = _b[dd]
	local diff_`ver'_dd = `r3_`ver'_dd'-`r1_`ver'_dd'
	
	qui suest r1_`ver' r3_`ver', vce(cluster geo_id)
	test [r1_`ver'_mean]dd = [r3_`ver'_mean]dd
	local diff_`ver'_dd_p = r(p)
	local diff_`ver'_dd_p : di %9.0g `diff_`ver'_dd_p'
	
	cap n qui reg `ver' dd i.geo_idnum i.yq [aweight=pop_FFS_2010] if cba_round==0, vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010] if e(sample)==1
		cap n outreg2 using $resdir\reg.dta, dta keep(dd) addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean), Panel C - Panel A, `diff_`ver'_dd', P value, `diff_`ver'_dd_p') `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA11_DIDCPAP_Adjacent") sheetreplace	
restore
}


* Robustness 
{
* Robustness 1: no non-CBA
local rp replace
foreach ver in price qX {
	cap n qui reghdfe lg_`ver' dd if cba_round!=0 [aweight=pop_FFS_2010], nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' if e(sample) [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}

* Robustness 2: no small non-CBA
preserve 
use ".\temp\ZIP_CBA_CBSA.dta", clear
keep geo_id insample 
duplicates drop 
tempfile temp 
save `temp'
restore
merge m:1 geo_id using `temp', keep(3) nogen
foreach ver in price qX {
	cap n qui reghdfe lg_`ver' dd if insample==1 [aweight=pop_FFS_2010], nocons absorb(geo_idnum yq) vce(cluster geo_id)
		qui sum `ver' if e(sample) [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean))
}

* Robustness 3: State FEs
preserve
use ".\temp\ZIP_CBA_CBSA.dta", clear
keep geo_id CBA state state_name cba_round
drop if state=="" | geo_id==""
drop if geo_id=="CBA15" & state_name=="ME" // CBA = Boston, MA-NH only 1 zipcode
duplicates drop
gen R1 = 1 if cba_round==1
gen R2 = 1 if cba_round==2
bysort state: egen R1X = max(R1)
bysort state: egen R2X = max(R2)
drop R1 R2
expand 28
bysort geo_id state: gen yq = yq(2008,4)+_n
format yq %tq
levelsof state, local(levels)
foreach l of local levels {
	gen state_`l' = (state=="`l'")
	gen state_`l'_post = 0
	replace state_`l'_post = 1 if state=="`l'" & R1X==1 & yq>=yq(2011,1)
	replace state_`l'_post = 1 if state=="`l'" & R2X==1 & yq>=yq(2013,3)
}
drop state_name
collapse (max) state_*, by(geo_id yq)
tempfile temp_state
save `temp_state'
restore
merge 1:1 geo_id yq using `temp_state', nogen

foreach ver in price qX {
	cap n qui reghdfe lg_`ver' dd state_*_post [aweight=pop_FFS_2010], nocons absorb(geo_idnum yq) vce(cluster geo_id)
	qui sum `ver' if e(sample) [aweight=pop_FFS_2010]
	cap n outreg2 using $resdir\reg.dta, dta keep(dd) addtext(Yr-Qtr FE, X, GEO FE, X, State FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean))
}


* Robustness 4: 28-quarter Enrollees
replace q2 = 0 if q2==.
gen q2X = q2/pop_FFS_full28Q*1000
gen lg_q2X = log(q2X)

replace n_bene2 = 0 if n_bene2==.
gen n_bene2X = n_bene2/pop_FFS_full28Q*1000
gen lg_n_bene2X = log(n_bene2X)

replace n_initiator2 = 0 if n_initiator2==.
gen Ninitiator2X = n_initiator2/pop_FFS_full28Q*1000
gen lg_Ninitiator2X = log(Ninitiator2X)	

replace total_amt2 = 0 if total_amt2==.
gen payment2X = total_amt2/pop_FFS_full28Q*1000
gen lg_payment2X = log(payment2X)

gen price2 = total_amt2/q2
gen lg_price2 = log(price2)

foreach ver in price2 q2X {
	cap n qui reghdfe lg_`ver' dd [aweight=pop_FFS_2010], nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean))
}

preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA18_DIDCPAP_robust") sheetreplace	
restore
}


*** Event Study Plot (CPAP) ***
{
use ".\temp\sample_CPAP_reg.dta", clear
gen CBP = .
replace CBP = yq(2011,1) if cba_round==1
replace CBP = yq(2013,3) if cba_round==2 & CBP==.
replace CBP = yq(2014,1) if cba_round==1 & CBP==.
gen relative_qtr = yq-CBP
tab relative_qtr, m
drop if relative_qtr>9 & relative_qtr!=.
forvalues i = -20/9 {
	local q = `i'+21
	gen Q`q' = (relative_qtr==`i')
}
egen geo_idnum = group(geo_id)
local qtr_list Q1 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11 Q12 Q13 Q14 Q15 Q16 Q17 Q18 Q19 /// Q20 as baseline
	 Q21 Q22 Q23 Q24 Q25 Q26 Q27 Q28 Q29 Q30 
local absorb_list "geo_id yq"

foreach ver in lg_price lg_qX lg_n_beneX lg_NinitiatorX lg_paymentX HHI_NPI N_NPI {
	if "`ver'"=="lg_price" {
		local t "Log Price"
		local v "Fig1"
	}
	if "`ver'"=="lg_qX" {
		local t "Log Quantity"
		local v "Fig1"
	}
	if "`ver'"=="lg_n_beneX" {
		local t "Log Beneficiaries"
		local v "Fig1"
	}
	if "`ver'"=="lg_NinitiatorX" {
		local t "Log New Beneficiaries"
		local v "Fig1"
	}
	if "`ver'"=="lg_paymentX" {
		local t "Log Spending"
		local v "Fig1"
	}
	if "`ver'"=="HHI_NPI" {
		local t "HHI"
		local v "FigA1"
	}
	if "`ver'"=="N_NPI" {
		local t "Number of supplier (w/ >1% market share)"
		local v "FigA1"
	}
	if "`ver'"=="top_share_payment" {
		local t "Top National Supplier Market Share"
		local v "FigA1"
	}

	qui reghdfe `ver' `qtr_list' [aweight=pop_FFS_2010], nocons absorb(`absorb_list') vce(cluster geo_id)
	outreg2 using $resdir\reg.dta, dta fmt(f) stat(coef se ci_low ci_high) noparen noaster replace
	
	preserve
	use $resdir\reg_dta.dta, clear
	gen i = 1 if v1!=""
	gen n = sum(i)
	bysort n: gen t = _n
	drop i
	reshape wide v1 v2, i(n) j(t)
	rename v11 var
	rename v21 coef
	rename v22 se
	rename v23 ci_low
	rename v24 ci_high
	keep var coef se ci_low ci_high
	gen relative_qtr = substr(var,2,.)
	keep relative_qtr coef se ci_low ci_high
	destring relative_qtr coef se ci_low ci_high, replace force
	replace relative_qtr = relative_qtr-21
	drop if relative_qtr==.
	set obs 28
	replace relative_qtr = -1 if relative_qtr==.
	foreach s in se coef ci_low ci_high {
		replace `s' = 0 if relative_qtr==-1
	}
	tostring relative_qtr, gen(qtr)
	label define relative_qtr_v -20 "-20" -16 "-16" -12 "-12" -8 "-8" -4 "-4" 0 "0" 4 "4" 8 "8", replace
	label value relative_qtr relative_qtr_v
	sort relative_qtr
	scatter ci_high ci_low coef relative_qtr, connect(direct direct direct) ///
		lc(black black black) ms(none none none) lp(dot dot solid) ///
		xline(0, lwidth(vthin) lp(dash)) xlabel(-20(4)8, valuelabel) legend(off) ///
		xtitle("Relative Qtr to CBP") title("`t'") graphregion(color(white))
	graph export $resdir\\`v'_`ver'_CPAP.png, as(png) replace
	graph export $resdir\\`v'_`ver'_CPAP.pdf, as(pdf) replace
	restore
}

}



*******************
*** Sleep Apnea ***
*******************

*** Sample construction ***
{
use ".\statadata\OSA_diag", clear
keep if year(date)>=2009
gen yq = qofd(date)
collapse (max) osa sa, by(bene_id yq)
save ".\temp\sample_OSA.dta", replace

gen diag_yq_OSA = yq if osa==1 
gen diag_yq_SA = yq if sa==1
collapse (min) diag_yq*, by(bene_id)
format diag_yq* %tq

* MBSF
merge 1:m bene_id using ".\temp\mbsf", keep(3) nogen
keep bene_id year geo_id cba_round enroll* dual* diag_yq*
reshape long enroll_Q dual_Q, i(bene_id year geo_id cba_round diag_yq*) j(qtr)
gen yq = yq(year,qtr)
format yq %tq

keep if enroll_Q==3 
gen dual = 1 if dual_Q==3
replace dual = 0 if dual_Q==0

merge 1:1 bene_id yq using ".\temp\sample_OSA.dta", keepusing(osa sa) keep(1 3) nogen

format diag_yq* %tq
egen bene_idnum = group(bene_id)
xtset bene_idnum yq

gen n_diag_OSA = 1 if diag_yq_OSA==yq
gen n_diag_acc_OSA = 1 if diag_yq_OSA<=yq
gen n_diag_tt_OSA = 1 if osa==1
gen n_diag_acc4_OSA = 1 if osa==1 | L.osa==1 | L2.osa==1 | L3.osa==1

gen n_diag_SA = 1 if diag_yq_SA==yq
gen n_diag_acc_SA = 1 if diag_yq_SA<=yq
gen n_diag_tt_SA = 1 if sa==1
gen n_diag_acc4_SA = 1 if sa==1 | L.sa==1 | L2.sa==1 | L3.sa==1

gen n_diag_SAnOSA = 1 if diag_yq_SA==yq & diag_yq_OSA>yq
gen n_diag_OSAaSA = 1 if diag_yq_OSA==yq & diag_yq_SA<yq
gen s_diag_OSAaSA = (n_diag_OSAaSA==1) if diag_yq_OSA==yq
gen n_diag_OSAbSA = 1 if diag_yq_OSA==yq & diag_yq_SA>=yq
save ".\temp\sample_OSA.dta", replace

use ".\temp\sample_CPAP.dta", clear
collapse (min) fst_CPAP_yq = yq, by(bene_id)
merge 1:m bene_id using ".\temp\sample_OSA.dta", keep(2 3) nogen
gen n_diag_OSAaCPAP = 1 if diag_yq_OSA==yq & fst_CPAP_yq<yq
gen s_diag_OSAaCPAP = (n_diag_OSAaCPAP==1) if diag_yq_OSA==yq
tab s_diag_OSAaCPAP
save ".\temp\sample_OSA.dta", replace

use ".\statadata\OSA_cost_sum.dta", clear
replace yq = qofd(yq)
format yq %tq
merge 1:1 bene_id yq using ".\temp\sample_OSA.dta", keep(2 3) nogen
foreach ver in payment payment_ipop payment_dme payment_d {
	replace `ver' = 0 if `ver'==.
	replace `ver' = . if `ver'<0
	gen `ver'_acc = `ver' if n_diag_acc_OSA==1
	gen `ver'_acc4 = `ver' if n_diag_acc4_OSA==1
}
gen payment_ab = payment-payment_d-payment_dme // payment AB without DME
replace payment_ab = . if payment_ab<0
gen payment_ab_acc = payment_ab if n_diag_acc_OSA==1
gen payment_ab_acc4 = payment_ab if n_diag_acc4_OSA==1
save ".\temp\sample_OSA.dta", replace

* regression sample
use ".\temp\sample_OSA.dta", clear
collapse (sum) n_diag* (mean) s_diag* payment*, by(yq geo_id cba_round) 
drop if yq==.
drop if geo_id==""
unique yq
unique geo_id
di 969*28
save ".\temp\sample_OSA_reg.dta", replace

keep geo_id cba_round
duplicates drop
expand 28
bysort geo_id: gen yq = yq(2008,4)+_n
format yq %tq
merge 1:1 geo_id yq using ".\temp\sample_OSA_reg.dta", nogen
merge 1:1 geo_id yq using ".\temp\population", keep(1 3) nogen

foreach v of varlist n_diag* {
	replace `v' = 0 if `v'==.
	gen `v'_X = `v'/pop_FFS*1000
	gen lg_`v'_X = log(`v')
}

foreach ver of varlist payment* {
	gen lg_`ver' = log(`ver')
}

gen dd = 0
replace dd = 1 if yq>=yq(2011,1) & yq<=yq(2013,4) & cba_round==1
replace dd = 1 if yq>=yq(2013,3) & yq<=yq(2016,2) & cba_round==2
replace dd = 1 if yq>=yq(2014,1) & yq<=yq(2016,4) & cba_round==1
save ".\temp\sample_OSA_reg.dta", replace
}


*** Regression: DID OSA ***
{
use ".\temp\sample_OSA_reg.dta", clear
egen geo_idnum = group(geo_id)

local rp replace
foreach ver in n_diag_OSA_X payment_ab_acc4 payment_dme_acc4 {
	cap n qui reghdfe lg_`ver' dd [aweight=pop_FFS_2010], nocons absorb(geo_id yq) vce(cluster geo_id)
	qui sum `ver' if e(sample) [aweight=pop_FFS_2010]
	cap n outreg2 using $resdir\reg.dta, dta addtext(Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA13_OSA") sheetreplace	
restore

* CBA YQ trend
gen CBA_yq = yq*(cba_round!=0)
local rp replace
foreach ver in n_diag_OSA_X payment_ab_acc4 payment_dme_acc4 {
	cap n qui reghdfe lg_`ver' dd CBA_yq [aweight=pop_FFS_2010], nocons absorb(geo_id yq) vce(cluster geo_id)
	qui sum `ver' if e(sample) [aweight=pop_FFS_2010]
	cap n outreg2 using $resdir\reg.dta, dta addtext(CBA time trend, X, Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("TableA13RR_OSA_CBAtrend") sheetreplace	
restore
}

* Main results with time trend
{
use ".\temp\sample_CPAP_reg.dta", clear
egen geo_idnum = group(geo_id)

gen CBA_yq = yq*(cba_round!=0)
local rp replace
foreach ver in price qX n_beneX NinitiatorX paymentX {
	cap n qui reghdfe lg_`ver' dd CBA_yq [aweight=pop_FFS_2010], nocons absorb(geo_id yq) vce(cluster geo_id)
		qui sum `ver' [aweight=pop_FFS_2010]
		cap n outreg2 using $resdir\reg.dta, dta addtext(CBA time trend, X, Yr-Qtr FE, X, GEO FE, X, SE cluster, geo, weights, pop_2010) addstat(Mean of Dep Var, r(mean)) `rp'
	local rp
}
preserve
use $resdir\reg_dta.dta, clear
export excel using $resdir\Table.xls, sheet("Table2RR_DIDCPAP_CBAtrend") sheetreplace	
restore
}

	
*** OSA/SA Event Study ***
{
use ".\temp\sample_OSA_reg.dta", clear
gen CBP = .
replace CBP = yq(2011,1) if cba_round==1
replace CBP = yq(2013,3) if cba_round==2 & CBP==.
replace CBP = yq(2014,1) if cba_round==1 & CBP==.

gen CBA_yq = yq*(cba_round!=0)
		
gen relative_qtr = yq-CBP
tab relative_qtr, m
drop if relative_qtr>9 & relative_qtr!=.
forvalues i = -20/9 {
	local q = `i'+21
	gen Q`q' = (relative_qtr==`i')
}
egen geo_idnum = group(geo_id)

local qtr_list Q3 Q4 Q5 Q6 Q7 Q8 Q9 Q10 Q11 Q12 Q13 Q14 Q15 Q16 Q17 Q18 Q19 /// Q20 as baseline
	 Q21 Q22 Q23 Q24 Q25 Q26 Q27 Q28 Q29 Q30 

foreach ver in OSA SA {

	if "`ver'"=="OSA" {
		local yver lg_n_diag_`ver'_X
		local ytt "Obstructive Sleep Apnea (OSA)"
	}
	if "`ver'"=="SA" {
		local yver lg_n_diag_`ver'_X
		local ytt "Sleep Apnea (SA)"
	}
	qui reghdfe `yver' `qtr_list' [aweight=pop_FFS_2010], nocons absorb(geo_id yq) vce(cluster geo_id)
		outreg2 using $resdir\reg.dta, dta fmt(f) stat(coef se ci_low ci_high) noparen noaster replace
	
	preserve
	use $resdir\reg_dta.dta, clear
	gen i = 1 if v1!=""
	gen n = sum(i)
	bysort n: gen t = _n
	drop i
	reshape wide v1 v2, i(n) j(t)
	rename v11 var
	forvalues v = 2/2 {
		rename v`v'1 coef`v'
		rename v`v'2 se`v'
		rename v`v'3 ci_low`v'
		rename v`v'4 ci_high`v'
	}
	gen relative_qtr = substr(var,2,.)
	keep relative_qtr coef* se* ci_low* ci_high*
	destring relative_qtr coef* se* ci_low* ci_high*, replace force
	replace relative_qtr = relative_qtr-21
	drop if relative_qtr==.
	set obs 28
	replace relative_qtr = -1 if relative_qtr==.
	foreach v of varlist coef* se* ci_low* ci_high* {
		replace `v' = 0 if `v'==.
	}
	tostring relative_qtr, gen(qtr)
	label define relative_qtr_v -20 "-20" -16 "-16" -12 "-12" -8 "-8" -4 "-4" 0 "0" 4 "4" 8 "8"
	label value relative_qtr relative_qtr_v
	sort relative_qtr
	scatter ci_high2 ci_low2 coef2 relative_qtr, connect(direct direct direct) ///
		lc(black black black) ms(none none none) lp(dot dot solid) ///
		xline(0, lwidth(vthin) lp(dash)) xlabel(-20(4)8, valuelabel) legend(off) ///
		xtitle("Relative Qtr to CBP") title("`ytt'") graphregion(color(white))
	graph export $resdir\FigA3_`ver'.png, as(png) replace
	graph export $resdir\FigA3_`ver'.pdf, as(pdf) replace
	restore
}

}


*** Dual vs. Non-Dual ***
{
use ".\temp\sample_OSA.dta", clear

collapse (sum) n_diag*, by(yq geo_id cba_round dual) 
drop if dual==.
drop if yq==.
drop if geo_id==""
unique yq
unique geo_id
di 969*28*2
tempfile temp
save `temp'
keep geo_id cba_round
duplicates drop
expand 28
bysort geo_id: gen yq = yq(2008,4)+_n
format yq %tq
expand 2
bysort geo_id yq: gen dual = _n-1
merge 1:1 geo_id yq dual using `temp', nogen
merge m:1 geo_id yq using ".\temp\population", keep(1 3) nogen

foreach v in OSA SA {
	replace n_diag_`v' = 0 if n_diag_`v'==.
	gen n_diagX_`v' = n_diag_`v'/pop_FFSdual*1000 if dual==1
	replace n_diagX_`v' = n_diag_`v'/pop_FFSnondual*1000 if dual==0
	gen lg_n_diagX_`v' = log(n_diagX_`v')
	replace n_diag_acc_`v' = 0 if n_diag_acc_`v'==.
	gen n_diag_accX_`v' = n_diag_acc_`v'/pop_FFSdual*1000 if dual==1
	replace n_diag_accX_`v' = n_diag_acc_`v'/pop_FFSnondual*1000 if dual==1
	gen lg_n_diag_accX_`v' = log(n_diag_accX_`v')
}

* Event Study Plot
egen geo_idnum = group(geo_id)
cap gen CBP = .
replace CBP = yq(2011,1) if cba_round==1
replace CBP = yq(2013,3) if cba_round==2 & CBP==.
replace CBP = yq(2014,1) if cba_round==1 & CBP==.
cap gen relative_qtr = yq-CBP
tab relative_qtr, m
drop if relative_qtr>9 & relative_qtr!=. 
forvalues s = -20/9 {
	local q = `s'+21
	cap gen Q_`q' = (relative_qtr==`s')
	cap gen Q1_`q' = (relative_qtr==`s' & dual==1)
	cap gen Q0_`q' = (relative_qtr==`s' & dual==0)
}

local qtr_list Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_21 Q1_22 Q1_23 Q1_24 Q1_25 Q1_26 Q1_27 Q1_28 Q1_29 Q1_30  /// Q20 as baseline
	Q0_2 Q0_3 Q0_4 Q0_5 Q0_6 Q0_7 Q0_8 Q0_9 Q0_10 Q0_11 Q0_12 Q0_13 Q0_14 Q0_15 Q0_16 Q0_17 Q0_18 Q0_19 Q0_21 Q0_22 Q0_23 Q0_24 Q0_25 Q0_26 Q0_27 Q0_28 Q0_29 Q0_30 // Q20 (relative_qtr==-1 as baseline)

foreach ver in OSA SA {
	
	if "`ver'"=="OSA" {
		local ytt "Obstructive Sleep Apnea (OSA)"
	}
	if "`ver'"=="SA" {
		local ytt "Sleep Apnea (SA)"
	}
	
	qui reghdfe lg_n_diagX_`ver' `qtr_list' [aweight=pop_FFS_2010], nocons absorb(geo_idnum#dual yq#dual) vce(cluster geo_id)
	outreg2 using $resdir\reg.dta, dta fmt(f) stat(coef se ci_low ci_high) noparen noaster replace
	
	preserve
	use $resdir\reg_dta.dta, clear
	gen i = 1 if v1!=""
	gen n = sum(i)
	bysort n: gen t = _n
	drop i
	reshape wide v1 v2, i(n) j(t)
	rename v11 var
	rename v21 coef
	rename v22 se
	rename v23 ci_low
	rename v24 ci_high
	gen relative_qtr = substr(var,4,.)
	gen group = substr(var,2,1)
	keep group relative_qtr coef se ci_low ci_high
	destring group relative_qtr coef se ci_low ci_high, replace force
	replace relative_qtr = relative_qtr-21
	drop if relative_qtr==.
	set obs 56
	replace relative_qtr = -1 if relative_qtr==.
	replace group = 1 in 55
	replace group = 0 in 56
	foreach s in se coef ci_low ci_high {
		replace `s' = 0 if relative_qtr==-1
	}
	tostring relative_qtr, gen(qtr)
	label define relative_qtr_v -20 "-20" -16 "-16" -12 "-12" -8 "-8" -4 "-4" 0 "0" 4 "4" 8 "8"
	label value relative_qtr relative_qtr_v
	sort group relative_qtr
	twoway (connect coef relative_qtr if group==1, ms(o) mc(gs10) lc(gs10)) ///
		(rcap ci_low ci_high relative_qtr if group==1, lc(gs10)) ///
		(connect coef relative_qtr if group==0, ms(t) mc(black) lc(black) lp(dash)) ///
		(rcap ci_low ci_high relative_qtr if group==0, lc(black)), ///
		xline(0, lwidth(vthin) lp(dash)) xlabel(-20(4)8, valuelabel) `ylabel' ///
		legend(order(1 3) lab(1 "Dual") lab(3 "Non-Dual")) ///
		xtitle("Relative Qtr to CBP") title("`ytt'") graphregion(color(white))
	graph export $resdir\\FigA6_`ver'_dual.png, as(png) replace
	graph export $resdir\\FigA6_`ver'_dual.pdf, as(pdf) replace
	restore
}

}














